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Edge reconstruction modifies the electronic properties of finite graphene samples. We formulate 
a low-energy theory of the reconstructed zigzag edge by deriving the modified boundary condition 
to the Dirac equation. If the unit cell size of the reconstructed edge is not a multiple of three with 
respect to the zigzag unit cell, valleys remain uncoupled and the edge reconstruction is accounted 
for by a single angular parameter i9. Dispersive edge states exist generically, unless |$| = tt/2. We 
compute -0 from a microscopic model for the "reczag" reconstruction (conversion of two hexagons 
into a pentagon-heptagon pair) and show that it can be measured via the local density of states. In 
a magnetic field there appear three distinct edge modes in the lowest Landau level, two of which 
are counterpropagating. 

PACS numbers: 73.22.Pr, 68.35.B-, 72.80.Vp, 73.21.Hb 



I. INTRODUCTION 

The bulk electronic properties of graphene^ are modi- 
fied by edge effects in a small sample. A prominent exam- 
ple is a narrow ribbon of graphene which, depending on 
the exact lattice termination, is either gapped (semicon- 
ducting) or metallicP Edge states may form a flat band 
which favors spin polarization)^ and may have appli- 
cations in spintronicsP Scanning tunneling microscopy 
(STM) has provided considerable experimental support 
for these predicted edge effectsPtH 

The honeycomb lattice of graphene can be termi- 
nated along different directions, with the zigzag and the 
armchair termination having the smallest unit cell (see 
Fig. [1^,) . Recent microscopic calculations have indicated 
that these edges are unstable against a reconstruction of 
the hexagonal lattice structure which increases the size 
of the unit cellP^HH] In particular, Koskinen et a/P^have 
shown that the lowest energy is reached for the zz(57) 
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FIG. 1: (Color online) (a) Zigzag and armchair edge of the 
honeycomb lattice of graphene. (b) zz(57), or reczag recon- 
struction of the zigzag edge. The translation vector T of the 
various edges is indicated, as well as the Bravais lattice vec- 
tors Ri and R2 of the honeycomb lattice (with two atoms A 
and B in the unit cell). 



reconstruction of the zigzag edge: two adjacent hexagons 
convert into a pentagon- heptagon pair (see Fig. [Ija) . The 
stability of this so-called reczag edge has been confirmed 
by a variety of theoretical calculations^^ and they have 
been observed by transmission electron microscopy! 17 * 18 ^ 
Electronic properties of the reczag edge (and related 
reconstructions) have been studied using the difference 
equations obtained from a tight-binding Hamiltonian on 
the terminated latticeF^HUl j n thi s paper we propose 
an alternati ve ap proach based on the Dirac differen- 
tial equation j 22 l 23 l with edge reconstruction accounted for 
through a boundary condition. 24 The two approaches are 
equivalent at low energies, when the wave length is large 
compared to the lattice constant. One advantage of the 
approach based on the Dirac equation is that it contains 
fewer independent parameters than the full tight-binding 
Hamiltonian. Another advantage is that the boundary 
conditions are strongly constrained by symmetry, pro- 
viding a simple criterion for the existence of edge states 
and the presence or absence of intervalley scattering. 

We show that a broad class of edge reconstructions 
can be described by a boundary condition governed by a 
single angular parameter §. These boundaries cause no 
intervalley scattering and support dispersive edge states 
for \&\ 7^ 7r/2. The i?-class of boundary conditions in- 
cludes any edge reconstruction having a unit cell that is 
m times the size of a zigzag unit cell, with m not divisible 
by three. Most importantly, the reczag edge (m = 2) be- 
longs to the $-class. The value of 1? can be computed from 
a microscopic model (and we will carry out this calcula- 
tion), but we also show how it can be directly measured 
by STM via the local density of states. 

The paper is organized as follows: In Sec. |Tl]we begin 
by discussing the general form of the boundary condition 
for reconstructed graphene edges and show how discrete 
symmetries can be used to reduce the number of free pa- 
rameters to one single parameter (the i9-class boundary 
condition). We then focus in Sec. |III| on the particular 
case of the reczag boundary and compute the numerical 
value of $ from a tight-binding model. Sees. |IV| and [V] 
are devoted to a calculation of the electronic structure of 
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graphene terminated by reczag edges without and with 



magnetic field, respectively. We conclude in Sec. VI The 
Appendices contain details of the calculations, as well as 
a discussion of the effects of next-nearest-neighbor hop- 
ping and edge potentials on the zigzag boundary condi- 
tion (which also belongs to the $-class, having m = 1). 



II. BOUNDARY CONDITION FOR 
RECONSTRUCTED EDGES 

A. Tight-binding and Dirac Hamiltonian 

We describe the electronic structure of graphene using 
the tight-binding Hamiltonian, 



H = Y,u 3 \i) 01, 



(2.1) 



with one orbital |z) per atom. In the bulk we restrict 
ourselves to uniform nearest-neighbor hopping with value 
t. Only close to the edge we allow for a reconstruction 
of the honeycomb lattice and variations in the hopping 
amplitudes iy. 

In the low-energy limit and sufficiently far from the 
boundary, excitations with energy e obey the Dirac equa- 
tion 



where the Hamiltonian 



H = v f t ® {tr • p) 
acts on a four-component spinor wave function 



(2.2) 



(2.3) 



* = *2, *3, #4) = (*A, -«*b, i^B, -*a)- ( 2 - 4 ) 

Here and ^ x denote the wave amplitude on the 
X G {A, B} sublattice in the valley K and K' respec- 
tively. The Fermi velocity is denoted by vf and p = 
(—ihd x , —ihdy) is the two-dimensional momentum oper- 
ator. The matrices Tj and <jj are the Pauli matrices in 
valley and sublattice space respectively (with unit matri- 
ces To and ctq)- 



The Dirac equation (2.2 1 has a sublattice (or "chiral") 
symmetry, 



fa 



)H{t z 



-H. 



(2.5) 



This symmetry implies that H i-> —H for \I/a | — ^ and 
— v£/ B . Physically, it expresses the fact that the 
nearest-neighbor hopping does not couple sites on the 
same sublattice. Chiral symmetry is preserved by lat- 
tice termination, but it is broken by edge reconstruction 
(which couples sites originating from the same sublat- 
tice). 



Boundary conditions for broken chiral 
symmetry 



The Dirac equation (2.2) must be supplemented by a 



boundary condition that also includes the effects of the 
edge reconstruction. 

In Ref. [24] it was shown that any valid current- 
conserving and time-reversally symmetric boundary con- 
dition for the Dirac equation has the form 



* = M*, M 



(n- 



n _L tib , (2.6) 



where ng is the unit vector in the x — y plane normal 
to the boundary, and u and n are three-dimensional unit 
vectors. If the edge makes an angle a with the x-axis, 
the boundary condition can be written more explicitly as 

\& = (y ■ t) ® (a z cos?? + (a x cos a + a y sin a) sin 

(2.7) 

with e (— 71-/2,71-/2]. 

Chiral symmetry requires that (t z ® <t z ) M(t z ® a z ) = 
M, which restricts the boundary condition ( 2.6 ) to zigzag 
[y = ±z, n = z) or armchair (y z = n z = 0) form. 
Since edge reconstruction breaks chiral symmetry, other 
boundary conditions are allowed. Still, we can reduce the 
three inde pend ent parameters of the general boundary 
condition ( |2.6| to a single parameter for a broad class of 
edge reconstructions, as we will now show. 

In the following we consider edges that are invariant 
under a lattice translation T = nR\ + mi? 2 , n,m G Z, 
where R x = (%/3a/2, -a/2) and R 2 = {Via/2, a/2) 
are the two Bravais lattice vectors of graphene. Fig. [l] 
shows the translation vector T for the example of the 
zigzag edge (n = — 1, m = 1), the armchair edge 
(n = 1, m = 1) and the reczag edge (n = —2, m = 2). 
Due to the translational symmetry the Bloch momentum 
k G [—7t/ |T| , 7r/ |T|] along the boundary is a conserved 
quantum number. A zone-folding argument, detailed in 
Appendix[X] shows that the two Dirac points of graphene 
project onto the same k if n = m mod 3 and different k 
otherwise. Conservation of k then implies that interval- 
ley scattering is forbidden unless n — m mod 3. 

These observations allow for some general statements: 
Any reconstruction of the armchair edge has a transla- 
tional vector T such that n = m mod 3, and hen ce al - 
lows for any three-parameter boundary condition (2.7 1. 
In contrast, any reconstruction of the zigzag edge has 
n = —to. Hence, if m is not divisible by 3, the boundary 
condition does not mix valleys. In this case v = ±z and 
the boundary condition (for a given edge orientation a) 
has the single-parameter form 

= ± t z ® (<J Z cos ■& + {a x cos a + <j y sin a) sin i9) \& , 
if n m mod 3. (2-8) 

The reczag boundary has a doubling of the unit cell 
with respect to zigzag (to = 2) and hence has boundary 
condition of the form (|2.8| . If however the unit cell is a 



tripled (or a multiple of a tripled) zigzag unit cell, the 
general boundary condition ( |2.7[ ) applies, i.e. valleys are 
typically mixed. An example of such an edge is the Z 2 n 
zigzag reconstruction discussed in Ref. 111! 

In the remainder of the paper we will focus on the 
reczag edge, since that has been predicted to be the most 
stable reconstruction.ESHUl However, we will give most of 
our results without specifying the angle <&, so that they 
apply to any edge with a boundary condition of the form 
(2.8). In order to emphasize this generality, we consider 



in Appendix [B] a zigzag edge where chiral symmetry is 
broken due to edge potentials or next-nearest-neighbor 
hopping, rather than due to edge reconstruction. 



III. BOUNDARY CONDITION FOR THE 
RECZAG EDGE 

A. Tight-binding model 



In order to obtain a value for the angle in Eq. ( 2.8 ) for 
the reczag edge, we employ a tight-binding parametriza- 
tion. We consider a reczag edge parallel to the y-axis 
(a — 90°), as shown in Fig. [2} The unreconstructed edge 
would have terminated with an atom of the B-sublatticc 
and we will therefore refer to the edge as the B-type 
reczag. (We give results for the A-type reczag at the end 
of the section.) The boundary condition for a B-type 
reczag edge along the y-axis reads 

* = -t z ® (era cos 1? + a y sin (3.1) 

We may write this boundary condition more explicitly in 



terms of the sublattice amplitudes (|2.4|) , 
with the definition 

T = tan(i?/2) 



$ 3 



(3.2) 
(3.3) 

(3.4) 



The reczag edge is translationally invariant over a dis- 
tance 2a, where a is the graphene lattice constant. Hence, 
wave functions in adjacent unit cells only differ by a phase 



/ 



r ,2ika 



with Bloch wave vector k. We allow for a vari- 



ation of the hopping amplitude due to the reconstruction, 
but assume for simplicity that the hopping amplitude on 
every hexagon remains given by the bulk value t. 

Numerical values for the modified hopping ampli- 
tudes from density functional theory (DFT) are in the 
literature^ (see Table [i]). An extended model for the 
reczag edge with more parameters has been studied in 
Ref. [21] We give results for the extended model in Ap- 
pendix [C] and show that there are no essential differences 
to the simpler model employed here. We also neglect the 
effects of hoppings beyond nearest-neighbor and edge po- 
tentials. These effects can all be accounted for by a mod- 
ification of the numerical value of i? (see Appendices |B| 
andlCl). 



<P*f 



L 1 \ / 




FIG. 2: (Color online) Nearest-neighbor tight-binding model 
of the reczag edge with identifiers for the hopping amplitudes 
(red) and the wave function amplitudes (blue). We take uni- 
form hopping amplitudes t away from the edge. The unit cell 
of the reczag edge is indicated in dark, the neighboring unit 
cells in light color. The wave functions in the neighboring 
unit cells are multiplied by a Bloch phase factor / = e 2lka . 
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0.99 


0.97 


1.5 


0.0753 


0.150 



TABLE I: DFT values for the hopping amplitudes t p in the 
tight-binding model for the reczag edge, from Ref. 1191 and 
the corresponding value of the boundary condition parameter 
T = tan(i?/2), calculated from Eq. \3.14\. 



Labeling wave function and hopping amplitudes as in- 
dicated in Fig. [2] we can write down the tight-binding 
equations, 



£fB 


= hipi + 


t((fA 4 


- f<PA) , 


(3.5a) 




= h<p2 + 


t (fA 4 


■<Pa), 


(3.5b) 




= h<pB -+ 


- ^2^2 4 


ft 3 tP4, 


(3.5c) 




= t 2 <pl + 


t3f3 + 


tltpB, 


(3.5d) 




= t 3 (p 2 + 






(3.5e) 


Eifii 


= hfi/f 


4 t 4 ip 3 




(3.5f) 



In the limit e — > it is now straightforward to find rela- 
tions for the wave functions on the first hexagons away 
from the reconstructed edge, 



<PA 



'■Pa 



ft\U 
(1 - f)t 



<Pb 



f>B 



m 



- ^2^4 
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£3 — ffati 

<Pb 
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tit 



21 1 



t 2 t4 



(3.6a) 
(3.6b) 



B. Boundary modes 

We proceed along the lines of Ref. [231 by separating 
the wave function ip into a part ^ that obeys the Dirac 
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equation, plus a boundary correction iphdy(r)- Since the 
valleys are not coupled, it is sufficient to consider a single 
valley at K = (0, K) = (0, ~4n/3a), 



Mr) = * B (r)e 



iK- 



bdyl 

+ ^dy(r). 



(3.7a) 
(3.7b) 



Taking further into account the translational symme- 
try along the y-direction we can write the wave function 
as 

Mr) = MjV Kv + $L y {jY k \ (3.8a) 
Mr) = MjV Kv + ^y{jV k \ (3.8b) 
K = K + n/a = -tt/3o. (3.8c) 

The index j numbers the unit cells transverse to the edge, 
with ifB, <£b corresponding to j = and ifA, if a to j = 1 
(see Fig. K2J). We denote by K the projection of the K- 
point into the doubled unit cell of the reczag edge. The 
Dirac modes thus have a periodicity given by the unper- 
turbed graphcnc lattice, whereas the boundary modes are 
governed by the periodicity of the recz ag re construction. 
Application of the boundary condition (3.1 1 on the Dirac 
modes specifies the angle $ € (— tt/2, 7r/2J from 



0a(O)/0b(O) =tan(0/2). 



(3.9) 



For the bulk of graphene away from the edge, the tight- 
binding equations take the form 

eip A (r) = t [Mr) + Mr - R%) + Mr - Ba)] , 

(3.10a) 

eMr) = t [Mr) + Mr + Ri) + M r + Ba)] • 

(3.10b) 



Inserting the decomposition (3.8) into Eq. (3.10) and ac 



counting for the fact that the Dirac and boundary modes 
have a different periodicity, we arrive in the limit e — > 
at 

MJ + 1) = Mj), tfbdyO' + 1) = ^Lyti), (3-lla) 

M3 + V = M3), ^ dy (i + l) = V3^ dy (i). (3.11b) 

In order for the wave function to be normalizable only 
non-growing contributions are allowed, so <^bdyC?) = 
for all j. The B-type reczag edge thus has a boundary 
mode on the A sublattice only. This boundary mode 
is a direct consequence of the unit cell doubling of the 
reconstructed edge. 

The boundary mode decays exponentially away from 
the edge, with a decay length of 3a/2. This is also the 
distance from the edge where the Dirac equation — which 
does not capture the boundary modes — is valid. Hence, 
the reczag edge can be faithfully treated within the Dirac 
approach, as there are deviations only within the first few 
unit cells away from the boundary. 



C. Boundary condition 

The wave amplitudes </?a,b and </?a,b near the reczag 
edge can be written in terms of the Dirac and boundary 
modes as 



VA=[^A(l)+^dy(l)] r l '\ 

£ A = [«Mi)-0 A dy (i)] r 3/4 , 

= <fe(°)> 

<pb = <mo)/- 1/2 . 

With this decomposition we find from Eq. (13.61) that 



(3.12a) 
(3.12b) 
(3.12c) 
(3.12d) 



0a(O) =T(fo(p), 
J" = tan(tf/2) 



2t(4 



(3.13) 
(3.14) 



The numerical values for J- and $ for the reczag edge are 
given in Table [TJ 

This concludes the derivation of the boundary condi- 
tion for the B-type reczag edge. For the A-type reczag, 
the role of the A and B sublattices is interchanged. We 
thus have the boundary conditions 



*i = iT *3 = -iF^4 



(3.15) 
(3.16) 



with the same value (3.14) of T ' ■ 

The zigzag boundary conditiorP corresponds to T = 
or T = oo. As one can see from Eq. (3.14), T vanishes 
if H 



<titl, so for these matched hopping amplitudes 
the doubling of the unit cell at the edge has no effect 
on the boundary condition. This explains why a zigzag- 
edge behavior was found in a tight-binding study of edge 
reconstruction for the special case that all hopping am- 
plitudes have their bulk values.^ 



IV. ELECTRONIC STATES 
A. Dirac solutions 

The knowledge of the boundary condition allows us to 
calculate electronic properties. In this section we con- 
sider zero magnetic field and then in the next section 
the effect of a magnetic field is included. Although we 
use the numerical values of the reczag edge obtained in 
the previous Section for plots and comparisons to tight- 
binding models, the analytical results we obtain are valid 
for arbitrary angles 

Since the reczag edge does not mix the valleys, it is pos- 
sible to consider the K and if'-points separately. From 
Eqs. (3.2) and (3.15) we see that, given a solution for a 



particular valley, substitution of T — > —l/J 7 gives a so- 
lution in the other valley. In what follows we focus our 
discussion on the K -point. 
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We consider either one or two reczag edges along the 
y-direction. The solution of the Dirac equation (2.3| at 
energy e has the form ty(x,y) — ip(x)e' lky with 



ip(x) = A 



hyp 



(k + iq) 



B 



(k - iq) 



-iqx 



(4.1) 



The wave vector k is real, q is real or imaginary, and the 
dispersion relation is e = ±/wp\/fc 2 + q 2 . The relative 
amplitudes A, B of the superposition have to be deter- 
mined by the boundary condition. 



a ) 0.1 



0.05 



-0.05 



-0.1 




b) 



0.2 



-2 -1.8 -1.6 

k [l/2a] 



o.i 



-0.1 



-0.2 




-2-1 1 2 
k [l/2a] 



B. Edge state dispersion 

To study the dispersion relation of the edge state we 
take a semi-infinite graphene sheet for x > 0, terminated 
with a B-type reczag edge at x = 0. 

We first focus on decaying solutions with an imaginary 
q = iz and energy |e| < |topfc|. These edge states are 
affected most prominently by the edge reconstruction. 
Keeping only the exponentially decaying part of (4.1) 
and substituting the boundary condition 
the equation 



find 



hvp(z — k) = Te. 
This only has a normalizable solution for 



l-T 2 

z = k — T = kcos-d > 0. 

1 + T 2 



(4.2) 



(4.3) 



The normalized edge state wave function then reads 



edge/ \ _ 



iP ease (x) 



cos 2 #/2 I 



with energy^ 2 ^ 
e(k) = 



-hvpksin'd, forfccosi9>0. 



(4.4) 



(4.5) 



The solution for the K' -valley is found by the replace- 
ment of J- — >• —1/J- in Eq. (4.2 1, yielding a solution with 
energy 



e(k) = hvp k sin ?9, for k cos $ < 0. 



(4.6) 



These edge states exist for any < tt/2. 

It is instructive to compare the recz ag edge state with 
the well-known zigzag counterpart J 2 ^ which corresponds 
to the limit i9 —> 0. In accord with the tight-binding 
calculationsj^the main difference between the two types 
of edge states is their energy dispersion: While the zigzag 
edge state features a dispersionless band e(k) — 0, the 
reczag edge state has a linear dispersion with velocity 
i>Fsin"#. This has implications for the density of states 
(see Sec. |IVC ). 

Furthermore, the zigzag edge state is exactly zero on 
one sublattice (the A sublattice for a B-type zigzag edge) , 



FIG. 3: (Color online) (a) Comparison between tight-binding 
(black circle) and Dirac equation results (solid lines) for the 
edge state dispersion near the A'-point. Results are shown 
for different values of the hopping amplitude £4 of the reczag 
edge. (All other hopping amplitudes are as in Table |I]) The 
continuum of bulk states in the Dirac cone is indicated in 
grey, (b) The data for £4 = 1.5t on a larger scale, showing 
both the K and K' points. 



whereas the reczag edge couples the two sublattices. The 
coupling is such that the two components of the wave 
function only differ by a constant factor, ipi(x) = Fip2{x) 
for all x, not only at the boundary. The wave function 
thus has the same decay length (fccosi?) -1 into the bulk 
on each sublattice!^ For $ —> tt/2 the decay length di- 
verges and the edge state disappears in the bulk. 



In Fig. [3^, we compare the edge state dispersion (4.5) 
with the results of the tight-binding model of the reczag 
edge. (The tight-binding results were calculated for a 
nanoribbon of width W = lOOOv^a, large enough that 
the opposite edges were essentially decoupled.) Results 
are shown for different values of J 7 , obtained by mod- 
ifying the value of £4 with respect to the DFT values 
in Table [TJ As expected, we find excellent agreement for 
small e, corresponding to A:- values close to the K or K' 
points. Away from these Dirac points, the two discon- 
nected edge states of the Dirac equation are connected 
by the tight-binding model, see Fig. [3b. 



C. Density of states 



To make contact with STM experiments, we calculate 
the local density of states (DOS) on sublattice j = A,B, 
given by 



Vfar) e„) (|(*nMr)| 2 + \(K)j(r) 



(4.7) 
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The sum runs over all eigenstates *& n , *S>' n in valley K, K' 
with energy e n . For the reczag edge state we find 



V^(e,x)=T 2 V c B a ^(e,x), 
V^(e,x) = M ; C °f ( ? u9( M )e- 2 " 



(4.8a) 



ttHvf I sin I 
with u = -e(hv F tan (4.8b) 

The function O(it) is the unit step function [O(tt) = 1 
for u > and zero otherwise] . The coefficients g s = g v = 
2 indicate the degeneracies due to the spin and valley 
degree of freedom. 

Integrating out the transverse coordinate x and sum- 
ming over both sublattices we find the total DOS per 
unit length of the edge, 



2? cdgc (e) 



2nhvp |sin$ 



0(u). 



(4.9) 



This result holds in the energy range |e| < (Hvf/o) \ sini?| 
(beyond which the Dirac equation breaks down) . Such a 
constant DOS was also found for the case that the edge 
state acquir es a dispersion due to next-nearest neigh- 
bor hopping! 27 ! 28 ! Compared to the zigzag case, where 
£> edge (e) ~ 6(e), the density of states is greatly reduced 
by the reconstruction, which may well prevent the ferro- 
magnetic instability of the zigzag edgeP 

In addition to the decaying edge state with imaginary 
q, there is a continuum of bulk states with real q. Then 
the term 



Hvf , , . s , . k- 

— (k + iq) =sgn(e)-^= 



iq 



Sgn(e)e^ (4.10) 



in Eq. (4.1| is a pure phase (sgn is the sign function). 



These bulk solutions are given by 
</> bulk (z) = C bull 



nik 



sin qx + sgn(e) F sin(gx -I- ip) \ 
i sgn(e) sin(ga; — ip) + iJ'sin qx I ' 

(4.11) 

with hv-pq — \e\ siny> > and normalization constant 

Cbuik = 7r _1 / 2 (l + T 2 + 2sgn(e) J-cos^) -172 . (4.12) 

The local DOS of the bulk states follows upon integra- 
tion, 



V^ lk (e,x) 



2irh 2 v 2 ? j 



dip \rf ulk (x) 



(4.13) 



For Jc 1 the integral can be evaluated analytically, 



vT\x,e) 



9s9v\£\ 

4:TTh 2 Vp 



l-J (0 + 2sgn( £ ).FJ 1 (0 



V 



+ T 2 [J (£)- J a (0] +0(^ 3 )), 
bulk f x s ) = gsgv l £ l ( i _ j 9 (fi) 



(4.14a) 



■ S gn(4F[j 3 (0 - Ji(0] 



(4.14b) 




ffl 

Q 0.2 



-0.1 0.1 

e x a/hvp 



FIG. 4: (Color online) Local density of states as a function of 
energy, at various distances from the reczag edge. The con- 
tributions from the A and B sublattices are shown separately 
in the two panels. The peak in the density of states evolves 
according to Eq. (4.151 (dashed line in lower panel). 



Away from the edge, 
j f approaches the ie-symmetric 



T>bulk 



with £ = 2x\e\/hv-p- 
2>B Ulk "> 3s5v|e|/2^ 2 «J 
DOS of an infinite graphene sheet. The boundary effects 
break this electron-hole symmetry, as a manifestation of 
the chiral symmetry breaking by the reczag boundary 
condition. 

Fig. [4] shows the full local DOS on each sublattice, 



The = V 



edge 



x>£ ulk with x e 



{A, B}. The edge state 
manifests itself as a peak in the local DOS on the B 
sublattice. The DOS on the A sublattice is much smaller 
near the edge (by a factor T 2 « 0.006). The peak energy 
£poak moves towards the Dirac point (the zero of energy) 
as the distance x from the edge is increased, according to 



^pcak 





2x' 



(4.15) 



for x > 3a/ 2, |#| <C 7r/2. (The Dirac approximation 
breaks down at smaller x, while for larger # the edge 
DOS no longer dominates over the bulk DOS.) We con- 
clude that STM experiments have direct access to the 
boundary condition angle through the dependence of 
the edge state peak on the distance from the edge. 



D. Nanoribbon 

So far we considered a semi-infinite graphene sheet 
with a single B-type reczag edge. Reczag nanoribbons 
(width W) will have a B-type reczag edge on one side (at 
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FIG. 5: (Color online) Comparison of the dispersion relations 
of modes in zigzag (a) and (b) reczag nanoribbons. The re- 
sults for both valleys are superimposed, by measuring k rela- 
tive to the if -point (black curves) and if '-point (blue curves), 
respectively. The dashed grey lines indicate the Dirac cone of 
graphene. 



x = 0) and an A-type reczag edge at the other side (at 
y = W). The spectrum now consists of a discrete set of 
transverse modes e n {k), governed by the transcendental 



equatio 



cos 2 i? (cos uj — cos 2 O) — sin 2 i9 cos u cos 2 
+ sin (sin — sin w sin 2d) = 0. 



(4.16) 



We defined lu 2 — ^W 2 [{e / hvp) 2 — k 2 ) andcosfi = Hvpk/e. 

Fig. [5] compares the mode dispersion of a zigzag 
nanoribbon 2 (i9 = 0) and a reczag nanoribbon. The 
prominent difference is the dispersion of the reczag edge 
mode: For kW cosd 3> 1 it is given (up to exponen- 
tially small corrections) by the results for a single edge, 
e(k) = ±hv-p |fc|sinz9, since then the wave functions on 
opposite edges decay rapidly and overlap only little. Both 
in the zigzag and reczag nanoribbon, the overlap of the 
edge states as k — > produces a larger and larger en- 
ergy splitting, until the edge states merge with the bulk 
bands. 

The bulk bands of the reczag nanoribbon have a slight 
offset towards towards negative energies (barely visible 
in Fig. [5]) , which breaks the electron-hole symmetry - 
again as a result of the breaking of chiral symmetry by 
the edge reconstruction. 



V. EFFECT OF A MAGNETIC FIELD 

A. Dirac solutions 

The presence of a uniform perpendicular magnetic field 
Bo is accounted for by the substitution p h-> p + eA, with 
— e the electron charge and A = B xy the vector poten- 
tial in the Landau gauge. The valleys remain uncoupled 
and translational invariance along the y-axis is preserved. 
The wave function ^>(x,y) = ?p(x)e lkv in a single valley 



thus satisfies the Dirac equation 
.V2 



in 
^2 



(5.1a) 
(5.1b) 
and Im = 



where E = el m /hvp, X = \[2 (x/Z m + kl n 
^h/eB is the magnetic length. 

The coupled first-order differential equations (5.1) de- 
couple into a second-order equation, 

d 2 x ^{x) = {\X 2 -\E 2 ±\)^{x), (5.2) 

where j — 1,2 and the plus sign holds for ijii while the mi- 
nus sign holds for ip2- Eq. (5.2) is solved by the parabolic 



cylinder function U{x,a) 1 determined up to normaliza- 
tion by^l 

d 2 x U = {\x 2 +a)U, ]imU(a,x) = 0. (5.3) 

x—i-oo 

The solution in a magnetic field takes the form 



E 
75 



1 — F 2 



ip 2 = iAU 



2 

1 + E 2 



BU 



,X 



iBU 



l-E 2 
1 + E 2 



-X 



(5.4a) 



-X 



(5.4b) 



where A and B are constants. 



B. Edge states and Landau levels 

Wc first consider a semi-infinite graphene sheet for x > 
0, terminated by a B-type reczag edge at x = 0. Only 
keeping the solutions that decay for x —± 00 in Eq. (5.4) 



and substituting the boundary condition (3.2), we obtain 



an implicit equation for the energy dispersion in the two 
valleys, 



E 

71 



u 



1+E- 1 



'2kL 



U 



' 1-E 2 



, V2kQ 



—T in valley K, 
l/T in valley K' . 



(5.5) 

The resulting dispersion is shown in Fig. [6] The main 
features can be understood from two principles: 

• The confining potential due to the magnetic field in 
Eq. (5.2) has its minimum at —kl 2 ^. Because of this, 



we find bulk-like Landau level solutions and hence 
flat bands for k <C with the bulk Landau level 
eneigy— e n = sgn{n){hv^/l m )^/2\n\, n G Z. For 
positive values of k, the center of the confining po- 
tential is moved beyond the edge of the sample, re- 
sulting in dispersive quantum Hall edge states with 
velocity vp (larger than the velocity v-p sin $ of the 
zero-field reczag edge states). 
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FIG. 6: (Color online) Energy dispersion at the B-type reczag 
edge in a magnetic field. The states in valley K and K' are 
shown in black and blue, respectively (with k measured rela- 
tive to the respective Dirac point). The zero-field edge state 
dispersion is included as dashed, red lines. 



• The magnetic field has little effect on the reczag 
edge states, if the edge state decay length is smaller 
than the magnetic length, |&;cos??| -1 -C l m . For 
this reason, we observe two bands in Fig. [6] that 
follow the reczag edge dispersion (shown as dashed 
lines) for large enough momenta. 



C. Triple edge mode in the lowest Landau level 

The interplay of the magnetic and zero- field edge states 
produces three distinct edge modes in the lowest Landau 
level (n = 0). These are labeled a,b,c in the top panel 
of Fig. [7| The unidirectional edge mode a in valley K is 
accompanied by a pair of counterpropagating edge modes 
in valley K' . These three modes have a distinct wave 
function profile, as shown in the lower panels of Fig. [7j 

For mode a in the if-valley, the bulk Landau level so- 
lution for k <C is nonzero on the B sublattice onlyPSl 
It moves closer to the edge with increasing k and even- 
tually becomes the reczag edge state, which is mostly 
localized on sublattice B, with a small 0(J- 2 ) contribu- 
tion on the A sublattice. In contrast, for modes b, c in 
the -RT'-valley, there are two solutions for every momen- 
tum: For k <C we find both the bulk Landau level 
solution (localized on sublattice A only) and the reczag 
edge state (localized mostly on sublattice B). Note that 
we find bands with a distinct bulk or edge character, in 
contrast to the zigzag edge where chiral symmetry forces 
always hybridized solutions.^ 

The tripling of the edge modes in the lowest Landau 
level does not change the value of the Hall conductance, 
since the contribution from the two counterpropagating 
modes cancels. But the valley polarization at the edge is 
changed. At a zigzag edge, the lowest Landau level edge 
modes are in the same valley for positive and negative 
energies, whereas they are in different valleys at an arm- 
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FIG. 7: (Color online) Nine lower panels: Probability den- 
sity profiles for three different values of k in the three distinct 
modes of the lowest Landau level, labeled a, b, c in the top 
panel. Mode a is in valley K, and the counter-propagating 
modes b, c are in valley K' . The colors distinguish the proba- 
bility densities on sublattice A (green) and B (red) . To allow 
a comparison of the profiles, the vertical axis in each graph 
has been rescaled. 



chair edgeF^l At the reczag edge both valleys are present 
for negative energy, with only a single valley for positive 
energy. 

D. Comparison with tight-binding model 

Fig. [8] shows a comparison between the band structure 
obtained from the Dirac equation and from the tight- 
binding model. (Similar tight-binding calculations are 
in Refs. [T9 21.) To be able to identify the contribu- 
tions from the two edges we took a wide nanoribbon, 
W = 8l m = 10lV3/2a, in which opposite edges are 
approximately decoupled. In this case the Dirac equa- 
tion results for the A-type reczag edge at x — W can 
be directly obtained from the results for a B-type reczag 
edge at x = by interchanging the valleys and replacing 
k-t-k-W/Q. 

The two calculations agree very well near the Dirac 
points. As in the zero-field case (Fig. |3Jd) the tight- 
binding model connects the edge states from the two val- 
leys K, K' , which are disconnected in the Dirac equation. 
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FIG. 8: (Color online) Comparison of the energy dispersion 
of a reczag nanoribbon in a magnetic field obtained from the 
tight-binding model (open circles) and from the Dirac equa- 
tion (red lines). For the lowest Landau level, the edge states 
localized at the x = and x — W boundary are highlighted in 
green and blue, respectively. The three lowest-Landau-level 
modes at the x = edge are labeled a,b,c . They appear 
displaced relative to Fig. [7J because there the momentum k 
is measured relative to the Dirac point of valley K, K' . 



VI. CONCLUSION 

In conclusion, we have derived the boundary condi- 
tion for the Dirac equation at reconstructed zigzag edges 
in graphene. The z9-class of boundary conditions (2.8) 



applies to reconstructions with a unit cell that is not a 
multiple of three times the zigzag unit cell. We have cal- 
culated the angular parameter $ for the zz(57) (reczag) 
reconstruction, which has been identified as the most sta- 
ble reconstruction. Most of our results are given for gen- 
eral < 7r/2, so they apply to other reconstructions in 
the $-class as well. 

The $-class reconstructions share two key properties: 
they do not cause intervalley scattering and they support 
edge states. Dispersive edge states were previously found 
for the reczag edgep' the zigzag edge with next-nearest 
neighbor hoppingpS and the zigzag edge with a bound- 
ary potential.^ Our analysis identifies an entire class of 
reconstructions with edge states, and gives analytic ex- 
pressions for the edge state dispersion in terms of a single 
parameter <&. 

The edge mode appears in the local density of states 
as a peak at energy e p0 ak- The dependence o f e pea k on 
the separation x from the edge, given by Eq. (4.15), al- 



lows a direct measurement of d by scanning tunneling 
microscopy. 

In a magnetic field there appears a tripling of the edge 
modes in the lowest Landau level. This could be ob- 
served in transport experiments, since two of three edge 
modes are counterpropagating and therefore susceptible 
to localization by disorder. With increasing disorder, the 



two-terminal conductance would then be reduced by a 
factor 1/3. 
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Appendix A: Condition for absence of valley mixing 
by edge reconstruction 

We explain the zone-folding argument used in Sec. 
|III A| to identify which periodicity of the edge recon- 
struction leaves the valleys uncoupled. It is similar to 
the zone-folding argument that distinguishes metallic and 
semiconducting carbon nanotubesP^ 

The projection of the i\-point along the direction of 
the edge is given by 



K ■ - — - = —(n—m) -, — r 

t r >\t\ 



and the projection of the iC'-point by 



K' 



1 27T 

-(771 — n — — ■ 
3 V '\T\ 



(Al) 



(A2) 



The projected K and If' -points correspond to the same 
momentum in the one-dimensional first Brillouin zone of 
the edge, if they differ by a multiple of a reciprocal lat- 
tice vector. This condition (K-K')-T/\T\ = l2n/\T\, 
I £ Z, is equivalent to the condition that n — m is di- 
visible by 3. Otherwise, if n ^ m mod 3, the if-points 
project to different momenta in the first Brillouin zone of 
the edge, and since these momenta are conserved due to 
translational symmetry, the valleys remain uncoupled. 



Appendix B: Boundary condition for modified 
zigzag edge 

Edge reconstruction is one modification of the zigzag 
edge that leads to a boundary condition of the single- 
In this Appendix we calculate the 



parameter form (2. 



value of the parameter i? for two alternative modifica- 
tions of the zigzag edge that break chiral symmetry: On- 
site potentials and next-nearest-neighbor hopping. Since 
most of our results for the reczag edge are given for ar- 
bitrary they can be applied to these edges as well — 
even though these modifications leave the lattice struc- 
ture unaffected. 

Consider a B-type zigzag edge with a nonzero potential 
Va, Vq on the outermost A and B atoms. (See Fig. [91.) 
Such on-site potentials could appear because the edge 



10 



a) 




; = o 1 



-1 
k [I /2a] 



FIG. 9: (Color online) (a) Schematic of the modified zigzag 
edge, with on-site potentials and hoppings labeled in red. (b) 
Comparison between the tight-binding edge state dispersion 
for the reczag edge (black circles), and the modified zigzag 
edge with Va = 0, ti = t, Vb = — Tt (blue squares). The 
Dirac equation has the same boundary condition at these two 
edges, leading to the same edge state dispersion near the Dirac 
point (red solid line). 



atoms see a different chemical environment than the bulk 
atoms. We also include a possible modification t\ of the 
hopping amplitude at the edge. The same model with 
Vb = —t 1 describes to leading order the effect of a next- 
nearest-neighbor hopping 

Since the unit cell is not changed by these modifica- 
tions, the boundary modes that appeared for the reczag 
edge are absent. Following the approach of Sec. |IIf| we 
find 



J" = tan(zV2) 



V A V B - t{ 



(Bl) 



This agrees with Rcfs. 28 34 for the special case Va = 
0, t\ — t. If next-nearest-neighbor hopping is the only 
modification, we set Vb = — f, Va. = 0, t% = t and arrive 
at 



J" = tan(tf/2) =t'/t. 



(B2) 



Fig. [9b shows a comparison of the edge state dispersion 



for the reczag edge from Sec. |III| and a zigzag edge with an 
edge potential such that the value of J- is the same. Both 



have the same boundary condition for the Dirac equation, 
and indeed we observe the same linearly dispersing edge 
state close to the Dirac point. 



Appendix C: Extended model for the reczag edge 

The tight-binding model for the reczag edge used in 
the main text is based on Ref. HH An extended model 
was studied in Ref. [STJ including also modifications of the 
hopping amplitudes in the first row of hexagons near the 
edge. From the general arguments of Sec.llljwe know that 




FIG. 10: (Color online) Schematic of the extended model for 
the reczag edge, with on-site energies and hoppings labeled in 
red. 



the form of the boundary condition remains the same, 
with a different numerical value for the parameter In 
this Appendix we calculate that value. 

The extended model of the reczag edge is shown in 
Fig. [10 In addition to the modified hopping amplitudes 
of Ref. 2H we also include (for additional generality) an 
on-site potential at the outermost edge atoms. Following 
the same procedure as in Sec. |III[ we obtain 



T = tan(tf/2) = T/Af, 



(CI) 



as the ratio of the coefficients 



T = tt\ [{t 2 (2t§ + 2t 5 *6 - 4) - 2{t\ + t 5 t 6 + t^Vx] (t\ ~ Vi) 



+t 2 3 {U(t\ - 2t 5 t 6 - 2tl) - 2(t 2 5 + tgta + tl)V 2 }\ , 



4 + {tl-v?){t 



v 2 2 ) + 4{t 2 t 4 



2ViV 2 )] 



(C2) 
(C3) 



Using the numerical values from Ref . [2U see Table pi} we find i? s» 0.0968 — within a factor of two from the 
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h/t 


h/t 


H/t 


U/t 




U/t 


J" 


1? 


0.94 


0.94 


1.06 


1.42 


1.04 


0.98 


0.0485 


0.0968 



value d 



0.150 following from the simpler model of Table 



TABLE II: Values of the hopping amplitudes in the extended 
tight-binding model, obtained from DFTp^ In this model, 
Vx = V2 = 0. The boundary condition parameter is calculated 



from Eq. (CI I 



1 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 

(2009) . 

2 L. Brey and H. A. Fertig, Phys. Rev. B 73, 235411 (2006). 

3 M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe, 
J. Phys. Soc. Jpn. 65, 1920 (1996). 

4 K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dressel- 
haus, Phys. Rev. B 54, 17954 (1996). 

5 Y.-W. Son, M. L. Cohen, and S. G. Louie, Nature 444, 
347 (2006). 

6 K. A. Ritter and J. W. Lyding, Nature Mat. 8, 235 (2009). 

7 Y. Niimi, T. Matsui, H. Kambara, K. Tagami, M. Tsukada, 
and H. Fukuyama, Phys. Rev. B 73, 085421 (2006). 

8 Y. Kobayashi, K. I. Fukui, T. Enoki, and K. Kusakabe, 
Phys. Rev. B 73, 125415 (2006). 

9 C. Tao, L. Jiao, O. V. Yazyev, Y.-C. Chen, J. Feng, X. 
Zhang, R. B. Capaz, J. M. Tour, A. Zettl, S. G. Louie, H. 
Dai, and M. F. Crommie, Nature Phys. 7, 616 (2011). 

10 P. Koskinen, S. Malola, and H. Hakkinen, Phys. Rev. Lett. 

101, 115502 (2008). 
1 T. Wassmann, A. P. Seitsonen, A. M. Saitta, M. Lazzeri, 

and F. Mauri, Phys. Rev. Lett. 101, 096402 (2008). 

12 B. Huang, M. Liu, N. Su, J. Wu, W. Duan, B.-L. Gu, and 
F. Liu, Phys. Rev. Lett. 102, 166404 (2009). 

13 J. Li, Z. Li, G. Zhou, Z. Liu, J. Wu, B.-L. Gu, J. Ihm, and 
W. Duan, Phys. Rev. B 82, 115410 (2010). 

14 G.-D. Lee, C. Z. Wang, E. Yoon, N.-M. Hwang, and K. M. 
Ho, Phys. Rev. B 81, 195419 (2010). 

15 C. K. Gan and D. J. Srolovitz, Phys. Rev. B 81, 125445 

(2010) . 

16 J. M. H. Kroes, M. A. Akhukov, J. H. Los, N. Pineau, and 
A. Fasolino, Phys. Rev. B 83, 165411 (2011). 

17 P. Koskinen, S. Malola, and H. Hakkinen, Phys. Rev. B 
80, 073401 (2009). 

18 C. Girit, J. C. Meyer, R. Erni, M. D. Rossell, C. 
Kisielowski, L. Yang, C.-H. Park, M. F. Crommie, M. L. 
Cohen, S. G. Louie, and A. Zettl, Science 323, 1705 (2009). 

19 P. Rakyta, A. Kormanyos, J. Cserti, and P. Koskinen, 
Phys. Rev. B 81, 115411 (2010). 

20 S. M.-M. Dubois, A. Lopez-Bezanilla, A. Cresti, F. Trio- 



zon, B. Biel, J.-C. Charlier, and S. Roche, ACS Nano 4, 
1971 (2010). 

21 J. N. B. Rodrigues, P. A. D. Goncalves, N. F. G. Rodrigues, 
R. M. Ribeiro, J. M. B. L. dos Santos, and N. M. R. Peres, 
larXiv: 1107.47791 

22 J. W. McChire, Phys. Rev. 104, 666 (1956). 

23 D. P. DiVincenzo and E. J. Mele, Phys. Rev. B 29, 1685 
(1984). 

24 A. R. Akhmerov and C. W. J. Beenakker, Phys. Rev. B 
77, 085423 (2008). 

25 V. A. Volkov and I. V. Zagorodnev, Low Temp. Phys. 35, 
2 (2009). 

26 Ref. I2U reaches a different conclusion, that the decay 
lengths into the bulk differ for the two sublattices in the 
case of a reczag edge. In their tight-binding approach the 
Dirac and boundary modes are not easily separated. While 
we find that the former have a sublattice independent de- 
cay length ~ 1/fc, the latter modes do introduce a sub- 
lattice dependence on the scale of the lattice constant, in 
accordance with Eq. ( |3.11 \. 

27 M. Wimmer, A. R. Akhmerov, and F. Guinea, Phys. Rev. 
B 82, 045409 (2010). 

28 J. Wurm, K. Richter, and I. Adagideli, Phys. Rev. B 84, 
075468 (2011). 

29 M. Abramowitz and I. A. Stegun, Handbook of Mathemat- 
ical Functions (Dover, New York, 1972). 

30 T. Ando, J. Phys. Soc. Jpn. 74, 777 (2005). 

31 L. Brey and H. A. Fertig, Phys. Rev. B 73, 195408 (2006). 

32 A. R. Akhmerov and C. W. J. Beenakker, Phys. Rev. Lett. 
98, 157003 (2007). 

33 K. Sasaki, S. Murakami, and R. Saito, Appl. Phys. Lett. 
88, 113110 (2006). 

34 S. Bhowmick and V. B. Shenoy, Phys. Rev. B 82, 155448 
(2010). 

35 R. Saito, G. Dresselhaus, and M. S. Dresselhaus, Physical 
Properties of Carbon Nanotubes (Imperial College, Lon- 
don, 2003). 

36 K. I. Sasaki, Y. Shimomura, Y. Takane, and K. Wak- 
abayashi, Phys. Rev. Lett. 102, 146806 (2009). 



